/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAWSWave
  \brief To calculate bound state radial wavefunction of a valence nucleon in
  a Woods-Saxon potential. Just for unit-test purposes.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022/03/07
  \date Last modified: 2022/03/07 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#ifndef _TAWSWave_h_
#define _TAWSWave_h_

#include <string>
#include <vector>
#include "TAODESolver.h"

using std::string;
using std::vector;

class TAWSWave : public TAODESolver{
public:
  /// inFile: take user input for potentials and single-particle state
  /// of the valence nucleon
  TAWSWave();
  virtual ~TAWSWave();

  double V(double r) const; ///< \retval V(r) = (V0+VSO+VC)*f2MuOverHbar2
  // potential for the unit-test method BoundExample
  double VT(double r) const; ///< \return total potential: V(r)+VL (VL is the centrifugal barrier)
  int Getl() const{ return l; }
  /// \retval Note that Rl=u*r; u = Rl/r. This method returns Rl, i.e. u*r, NOT u
  double Rl(double r); ///< implement an interpolation
  double E() const{ return fE; }

  /// the ODE set: du/dr=u1; du1/dr=-(lambda+v)*u; du2/dr=0,  where E=hbar^2/(2mu)*lambda ///
  /// the two point border condition: u(0)=0; du(0)/dr=1; u(\infty)=0; ///
  /// the ODE set itself: dyi/dx = f_i(x,y1,..yN)
  virtual void derivs(double x, const double *y, double *dydx) override;
  /// calculates the n-vector y[0..n-1] (satisfying the starting boundary conditions, of course)
  /// given the freely specifiable variables of v[0..n2-1] at the initial point x1
  // solve the radial wavefunction //
  void Bound(); ///< for solving the objective radial equation
  bool HasSolved(){ return fSolved; } ///< if Bound() has been called or not

private:
  bool fSolved; ///< if Bound() has been called or not

  /// single-particle state for the valence nucleon
  int n, l;
  double j;

  double *fxx, *fyy2; ///< the concatenated x and y^2
  int fcnt; // the length of the solution fyy
  double fE; ///< the solved eigen energy in MeV
};

#endif
